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ABSTRACT 

We explore a Boltzmann scheme for studying the evolution of compact binary 
populations of globular clusters. We include processes of compact-binary forma- 
tion by tidal capture and exchange encounters, binary destruction by dissociation 
and other mechanisms, and binary hardening by encounters, gravitational radi- 
ation and magnetic braking, as also the orbital evolution during mass transfer, 
following Roche lobe contact. For the encounter processes which are stochastic 
in nature, we study the probabilistic, continuous limit in this introductory work, 
deferring the specific handling of the stochastic terms to the next step. We focus 
on the evolution of (a) the number of X-ray sources Nxb in globular clusters, 
and (b) the orbital-period distribution of the X-ray binaries, as a result of the 
above processes. We investigate the dependence of Nxb on two essential cluster 
properties, namely, the star-star and star-binary encounter-rate parameters V 
and 7, which we call Verbunt parameters. We compare our model results with 
observation, showing that the model values of N X b and their expected scaling 
with the Verbunt parameters are in good agreement with results from recent 
X-ray observations of Galactic globular clusters, encouraging us to build more 
detailed models. 

Subject headings: globular clusters: general — binaries: close — X-rays: binaries 
— methods: numerical — stellar dynamics — scattering 

1. Introduction 

In this era of high-resolution X-ray observations with Chandra and XMM-Newton, stud- 
of compact binaries in globular clusters have reached an unprecedented level of richness 
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and detail. The numbers of compact X-ray binaries detected in Galactic globular clus- 
ters with high central densities are now becoming large enough that diagnostic correlations 
with essential cluste r parameters, such as the two-body encounter rate T, can be performed 
( IPooley et al. 1120031 ) at a high level of statistical significance. The results of such observa- 
tional studies are naturally to be compared with those obtained from theoretical modeling 
of binary dynamics in globular cl usters, which has had a long history, from the pioneering 
semi-analytic work of the 1970s (jHeggie II1975T). to the more detailed numerical scattering 



experiments of the 1980s (I Hut & BahcaJl 111983). leading to the wealth of detailed numeri 



cal work of the early- to mid-1990s (IMakino fc Aarseth Ill992l ; I Heggie &: Hut 1120031 ) using a 
variety of techniques including Fokker-Planck and Monte Carlo approaches, as also N-body 
simulations, and finally to the extensive N-body simula tions in the latter half of the 1990 s 
using special-purpose computers with ultrahigh speeds (IMakino fc Taiji 1 1 19981 : 1 Hut 1120011 ). 



The range of problems studied by the above modeling has also been extensive. From 
the study and classification of individual scattering events to th e construction of com 



prehensive fitting formulas for the cross-sections of such events (I Hut fc Bahcall 1 11983 



Heggie. Hut fc McMillan Ill996l ). from the development of Fokker-Planck codes to the use 



of M o nte Carlo methods for following b inary distributions in globular clusters (IGao et.al. 



1991 



Hut. McMillan &: Romani 1119921 ). and from tracking the fate of a relatively modest 
population of test binaries against a fixed stellar background to being able to tackle sim- 
ilar projec ts for much larger binary populat ions with the aid of the above special-purpose 



machines ( I Hut et.al. I Il992l ; IMakino I Il996l ) , efforts along various lines of approach have 
shed light on the overall phenomenon of binary dynamics and evolution in globular clus- 
ters from various angles. For example, evolutions of the distributions of both external and 
internal binding energies of the binaries under stellar encounters have been studied by sev- 
eral authors, the emphasis usually being on the former, and final results on the external 
binding energy being expressed almost universally in ter ms of their radial positions r in 



side the cluster, which provides an equivalent description ( 1 Hut. McMillan fc Romani 111992 



Sigurdsson &: Phinney 



19931 . 119951 ). 



In this series of papers, we introduce an alternative method of studying the evolution 
of compact-binary populations in globular clusters, wherein we use a Boltzmann description 
to follow the time-evolution of such populations, subject to both (a) those processes which 
determine compact-binary evolution in isolation {i.e., outside globular clusters, or, in the 
"field" of the host galaxy, so to speak), e.g., angular momentum loss by gravitational radia- 
tion and magnetic braking, as also orbital evolution due to mass transfer, and, (b) those pro- 



(Heggie 


1975; 


Shull 


1979: 



20061 ). binary formation through tidal capture and exchange processes, and binary destruc- 
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tion. We treat all of the above processes simultaneously through a Boltzmann formalism, 
the aim being to see their combined effect on the compact-binary population as a whole, 
in particular on the evolution of (a) the total number of X-ray binaries as the formation 
and destruction processes continue to operate, and, (b) the orbital-period distribution of the 
population. We stress at the outset that ours is not a Fokker-Planck description but the 
original Boltzmann one, which in principle is capable of handling both the combined small 
effects of a large number of frequent, weak, distant encounters and the individual large ef- 
fects of a small number of rare, strong, close encounters. In our approach, both of the above 
two types of effects are taken into account through cross- sect ions for the relevant processes, 
as determined from ex tensive previous work on numerical experiments with two-body and 



three-body encounters (I Heggie. Hut &: McMillan Ill996l ; l Portegies Zwart et.al. Ill997bl ). As 



these processes are inherently stochastic, a natural question that arises is how they are to be 
handled simultaneously with those which govern the fate of isolated compact binaries, and 
which are inherently continuous. It is essential to appreciate the importance of this ques- 
tion, since a simultaneous action of the above continuous and stochastic processes is precisely 
what operates on binaries in globular clusters, and so produces the observed properties of 
compact-binary populations in them. 

Our answer to the above question is a step-by-step one, as follows. As the first step, 
in this first paper of the series (henceforth Paper I), we explore the continuous limit of 
the above stochastic processes, wherein the probability or cross-section of a particular such 
process happening with a given set of input and output variables is treated as a continuous 
function of these variables. This is, of course, a simplification, but it serves as a clarification of 
the average, long-term trends expected in the evolution of the binary population. In the next 
step, in the second paper of the series (henceforth Paper II), we treat the stochastic processes 
as stochastic terms in the Boltzmann equation with cross- sect ions as given in Paper I, with 
the aid of relatively recently-developed methods for solving stochastic partial differential 
equations. The resulting evolutionary trends show stochastic behavior, as expected, with 
fluctuations that vary from one particular "realization" of the essential processes to another. 
However, the average trends follow the continuous limit computed in Paper I, which is as 
expected, and which shows the relevance of extracting this limit. 

In Papers I and II, we model the stellar background provided by the globular clus- 
ter as a fixed background with given properties, as has been widely done in previous works 



(I Hut. McMillan fc Romani II 1992 ; I Portegies Zwart et.al. Ill997bl ; ISigurdsson & Phinnev 111993 



19951 ): this amounts to neglecting the back reaction of binary evolution on the background, 
which is reasonable if the main aim is an investigation of essential features of binary evolution, 
as was the case in the above previous works, as also in this work. However, the globular- 
cluster background does evolve slowly, passes through the core-collapse phase and possible 
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gravothermal oscillations (jSugimoto fc Bettwieser 1 1 1 9 3 3l : iGao et.al. Ill99ll ). so that it would 



be interesting to be able to follow the effects of these on the evolution of the compact-binary 
population. We do this in the third paper of the series (henceforth Paper III), wherein 
we adopt previous results on time-evolution of globular-cluster properties, and study their 
effects on the evolution of compact-binary populations, again under the approximation of 
neglecting the back reaction of binary evolution on the globular-cluster background, as above 
and as appropriate for a first look. 

In our study, we focus primarily on two aspects of the compact-binary populations of 
globular clusters. First, we study how the total number Nxb of X-ray binaries (henceforth 
XBs, which are mass-transferring compact binaries where the donor is a low-mass "normal" 
star, and the accretor is a degenerate star — a neutron star or a heavy white dwarf) in 
a cluster evolves as the stellar encounter processes proceed. Second, we also follow the 
evolution of the orbital-period (P) distribution of the pre- X-ray binaries (henceforth PXBs; 
also see below) and XBs, (or, equivalently, the distribution of their orbital radii a) within 
the framework of our model. However, we have adopted here only a very simple model of 
orbital evolution of individual binaries in order to assess the feasibility of our basic approach 
to globular-cluster environments, as detailed later. Consequently, while the P-distribution 
found by us may be roughly applicable to cataclysmic variables (CVs) with white-dwarf 
accretors, it cannot be compared at this stage to that of low-mass X-ray binaries (LMXBs) 
with neutron-star accretors, without including the essential stellar evolutionary processes 
that occur during the PXB and XB phase. Thus, we record our computed P-distribution 
here only as a preliminary indication of the results that emerge naturally from this line of 
study at this stage, to be improved upon later. 

The basic motivation for our study comes from recent advances in X-ray observations 
of globular clusters, as mentioned above: with sufficient numbers of X-ray binaries detected 
in globular clusters, an understanding of how N X b is influenced by essential globular-cluster 
parameters is becoming a central question. With the above goal in mind, we therefore 
explicitly follow the evolution of binaries only in internal binding energy (or binary period, 
or binary separation, which are equivalent descriptions if the stellar masses are known) and 
time, but not of their external binding energy (or position inside the globular cluster; see 
above). We emphasize that we do not neglect changes in the latter in any way, as they are 
automatically taken care of in the detailed dynamics of encounters which are represented by 
the relevant cross-sections mentioned above and elaborated on in the following sections. It is 
only that we do not keep an explicit account of them, as we do not need them for our purposes. 
In other words, we consider a bivariate binary distribution function n(E in ,t), which may be 
looked upon as the integral of the distribution p(E ex , E in , t) over all admissible values of E KX , 



or equivalently over all positions r inside the globular cluster (1 Hut. McMillan fc Romani 
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19921 ; ISigurdsson fc Phinney I Il993l . Il995l ). We also emphasize that, by doing so, we do 



not implicitly assume any particu lar correlation, nor a lack thereof, between E in and E e 



(I Hut. McMillan fc Romani 1119921 ): whatever correlations result from the dynamics of the 
encounters will be automatically displayed if we follow the evolution in E ex or r, which is 
not of interest to us in this particular study. 

Our first results from the above evolutionary scheme show that the total number Nxb 
of XBs expected in a globular cluster scales in a characteristic way with well-known globular 
cluster parameters V and 7 (which we call Verbunt parameters: see Sec. 12.11) whose qualita- 



tive n ature is rather similar to that found in our earlier "toy" model (1 Banerjee fc Ghosh 



20061 ). although some details are different. Basically, Nxb scales with Y — a measure of 
the dynamical formation rate of compact binaries, and, at a given T, Nxb decreases with 
increasing 7 at large values of 7 — a measure of the rate of destruction of these binaries 
by dynamical processes. These expected theoretical trends with the Verbunt parameters 
compare very well with the observed trends in recent data, encouraging us to construct more 
detailed evolutionary schemes. 

In Sec. m we detail our model of compact binary evolution in globular clusters, describ- 
ing, in turn, our handling of globular clusters, binary formation, destruction, and harden- 
ing processes, our Boltzmann scheme for handling population-evolution, and our numerical 
method. In Sec. [31 we give our model results on (a) the expected number of X-ray bina- 
ries in globular clusters as a function of their Verbunt parameters, and (b) the evolution 
of compact-binary period distribution. In Sec. HI we compare these model results with the 
current observational situation. Finally, we collect our conclusions and discuss future possi- 
bilities in Sec. [51 



2. Model of Compact Binary Evolution in Globular Clusters 

We consider a binary population described by a number distribution n(a,t), where a is 
the binary separation, interacting with a fixed background of stars representing the core of a 
globular cluster of stellar density p and core radius r c . We now describe various ingredients 
of our model and the evolutionary scheme. 



2.1. Globular clusters 



Globular cluster cores are described by an average stellar density p, a velocity dispersion 
v c , and a core radius r c . In this work, we consider star-star and star-binary encounters 
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of various kinds, but neglect binary-binary encounters. F or characterizing the former two 
processes, two encounter rates are defined and used widely (I Verbunt 1120021 . 120061 ). The first 
is the two-body stellar encounter rate T, which scales with p 2 r^./v c , and occurs naturally in 
the rates of two-body processes like tidal capture, stellar collisions and merger. In fact, we 
can define it as 

r es ^ oc p^rl (1) 

V c 

for our purposes here. Note that the last scaling in the above equation holds only for 
virialized cores, where the scaling v c oc p l ^ 2 r c can be applied. In this work, we shall use this 
assumption where necessary, but with the caveat that some observed globular clusters have 
clearly not virialized yet. 

The second is a measure of the rate of encounter between binaries and single stars in the 
cluster, the rate normally used being the encounter rate 7 of a single binary with the stellar 
background, with the understanding that the total rate of binary-sin gle star encounter in 



the clu ster will be oc 717. We can define 7 for our purposes as we did in I Banerjee fc Ghosh 
(120061 ) . namely, 

(2) 



7 = — oc p l l 2 r~ 

v r 



where the last scaling holds, again, only for virialized cores. 

The importance of the above cluster parameters Y and 7 in this context has been 
extensively discussed by Verbunt ( I Verbunt 1 120021 . 120061 ). and we shall call them Verbunt 
parameters here. Note that, for virialized cores, we can invert Eqs. (CD) and (T5]) to obtain the 
scaling of the core density and radius with the Verbunt parameters as: 



p oc 



p2/5 7 4/5 



r c oc r 1/5 7~ 3/5 



(3) 



It is most instructive to display the observed globular clusters in the T — 7 plane, which 
we dc0 in Fig. [U The point that immediately strikes one in the figure is that the observed 
globular clusters seem to occur in a preferred, diagonal, "allowed" band in the T — 7 plane, 
along which there is a strong, positive correlation between the two parameters. We shall 
return to the significance of this elsewhere. 

In Fig. (U we also overplot the positions of those clusters in which significant numbers 
of X-ray sources have been detected, color-coding them according to the number of X-ray 
sources in each of them, as indicated. It is clear that these clusters are all in the upper parts 



1 Alternatively, the display can be in the p — r c plane, as in Verbunt's original work. We find the cluster 
dynamics more transparent when shown directly in terms of the Verbunt parameters. 
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of the above "allowed" band, which is entirely consistent with the widely-accepted modern 
idea that the dominant mechanisms for forming these compact XBs in globular clusters are 
dynamical, e.g., tidal capture, exchange encounters, and so on, since such mechanisms occur 
more efficiently at higher values of the Verbunt parameters T and 7, corresponding to higher 
stellar densities in the cluster core. Note that the probability of destruction of binaries by 
dynamical processes also increases with increasing 7, as we shall see below, so that, at first 
sight, we might have expected the highest incidence of XBs in those clusters which have high 
T and low 7. However, since T and 7 are strongly correlated positively, as above, we cannot 
have arbitrarily high T and low 7 for the same cluster. In reality, the highest number of 
XBs seem to occur, as Fig. [1] shows, in those clusters which have the highest values of V and 
high, but not the highest, values of 7. We return to this point later in the paper, where we 
present our theoretical expectations for the scaling of the number of binary X-ray sources 
with the Verbunt parameters T and 7 on the basis of the evolutionary scheme explored here. 

In modeling the globular cluster core as a static background in this work, we assume 
that, initially, a fraction kb of the stars is in primordial binaries, and that a fraction kx of 
the stellar population is compact, degenerate stars with the canonical mass nix — 1.4M Q 
(representing neutron stars and heavy white dwarfs). The rest of the stellar background 
(including the primordial binaries) is taken to consist of low-mass stars of the canonical 
mass rrif = 0.6 M(7), which is a reasonable estim ate of the mean stellar mass of a mass- 



seggregated core (I Portegies Zwart et.al. Ill997al ). Naturally, the compact binaries formed 



from these ingredients consist of a degenerate star of mass nix — 1.4M , and a low-mass 
companion of mass m c = rrif = O.6M . While this is clearly an oversimplification which 
must be improved upon in subsequent work, it appears to be adequate for a first look, which 
is our purpose here. 



2.2. A Boltzmann evolutionary scheme 

We explore in this work a Boltzmann evolutionary scheme, wherein the evolution of the 
number n(a, t) of binaries per unit interval in the binary separation a (we choose to work 
here with a; equivalent descriptions in terms of the binary period P or the internal binding 
energy [see Sec. [1] E in are possible, of course) is described by 

^ = R(a)-nD(a). (4) 

Here, Dn(a,t)/Dt = dn/dt + (dn/da)(da/dt) is the total derivative of bivariate n(a,t): as 
explained in Sec. [T], this n(a, t) is the result of an integration of a general, multivariate binary 
distribution over the variables we do not follow explicitly in this study, e.g., the external 



8 



binding energy or, equivalently, the position of the binary inside the globular cluster. Further, 
R(a) is the total rate of binary formation per unit interval in a due to the various processes 
detailed below, and D(a) is the total rate of binary destruction per binary per unit interval in 
a due to various processes, also detailed below. As our model stellar background representing 
the cluster core is taken as static for Papers I and II, the Verbunt parameters T and 7 are 
time-independent, so that the formation and destruction rates R and D only depend on a 
and the stellar masses. 

The above evolution equation can be re-written in the usual Boltzmann form 



where f(a) = da/dt represents the total rate of shrinkage or hardening of binaries (i.e., 
da/dt < 0) due to several effects, which we introduced in Sec. [Q and which we elaborate 
on below. In the absence of all processes of formation and destruction, R(a) = = D(a), 
Eq. ((SJ) becomes the usual collisionless Boltzmann equation 



representing a movement or "current" of binaries from larger to smaller values of a due to 
hardening. Equation (J6j) as akin to a wave equation with a formal "phase velocity" f(a) 
of propagation. This analogy often proves useful for solving many problems, even with the 
more complicated formation and destruction terms present in Eq. (jSJ). Note that, when 
f(a) is constant (or roughly so, which can happen under certain circumstances, as we shall 
see later), the elementary wave-equation analogy is quite exact, and solutions of the form 
n(a — f t) should apply. We shall explore this point elsewhere. 

Note further that the Boltzmann scheme outlined above does not have an explicit in- 
clusion of the escape of those binaries from the globular cluster which receive a sufficiently 
large "kick" . In principle, we can include this by suitably generalizing the above destruction 
term D(a). However, in this introductory study, this did not appear crucial, as the main 
population affected by this process is that of primordial binaries, whereas our main concern 
here is with dynamically-formed compact binaries. The latter are, generally speaking, al- 
ready so hard at formation that this process is much less effective in ejecting them from the 
cluster. Accordingly, we neglect this process here. 




(5) 




(6) 



2.3. Binary hardening processes 



In all of the dynamical encounter processes considered in this work, viz., collisional 
hardening (described in this subsection), and dynamical formation and destruction processes 
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(described in the next subsections), we shall assume the orbits to be circular, i.e., neglect 
their eccentricity. This is, again, a simplification used for a first look. However, it is well- 
known from extensive numerical simulati ons that a large majority of th e binaries formed 
by tidal capture are circular or nearly so (I Portegies Zwart et.al. Ill997bl ). due to the rapid 
circularization which follows capture. Since our main concern here is with dynamically- 
formed binaries, this approximation may well be a reasonable one for describing overall 
evolutionary properties of such binary populations. 



2.3.1. Hardening in pre-X-ray binary (PXB) phase 



As explained in detail in I Banerjee fc Ghosh I (120061 ). referred to henceforth as BG06 



the processes that harden binaries are of two types, viz., (a) those which operate in isolated 
binaries, and are therefore always operational, and (b) those which operate only when the bi- 
nary in a globular cluster. In the former category are the processes of gravitational radiation 
and magnetic braking, and in the latter category is that of collisional hardening. As dis- 
cussed in detail in BG06, collisional hardening, which increases with increasing a, dominates 
at larger orbital radii, while gravitational radiation and magnetic braking, which increase 
steeply with decreasing a, dominate at smaller orbital radii. It is these processes that harden 
a compact binary from its pre-X-ray binary (PXB) phase, during which its orbit is still not 
narrow enough for the companion (mass donor) star to come into Roche lobe contact, to 
the state where this Roche lobe contact does occur, at which point the companion starts 
transferring mass to the degenerate star, and the system turns on as an X-ray binary (XB) 
- either a CV or a LMXB, depending on the nature of the degenerate accretor. 

Consider gravitational radiation first. The relative angular momentum loss rate due to 
this process is: 

(j\ 32G 3 
-J \ = -a GW a~ 4 , a GW = —£-m c m x (m c + m x )- (7) 
/ GW 

Here, as before, mx is the mass in solar units of the degenerate primary (neutron star or 
white dwarf) which emits X-rays when accretion on it occurs during the mass-transfer phase 
of the compact binary, m c is the mass of its low-mass companion in solar units, and the unit 
of the binary orbital radius a is the solar radius. We shall use these units throughout the 
work. 



Now consider magnetic braking. The pioneering Verbunt-Zwaan (I Verbunt &: Zwaan 



198ll ) prescription for this process has been reassessed and partly revised in recent years, 



in view of further observational evidence on short-period binaries available now (for fur- 
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ther details, see discussions in BG06 and references therein), and modern prescriptions are 



suggested in I van der Sluys et.al. I (120051 ) . From these, we have chosen for this work the 
following one which preserves the original Verbunt-Zwaan scaling, but advocates an overall 
reduction in the strength of the magnetic braking process: 

M 3 

3mb(o) = | ^ | = -a MB a 5 , a MB = 9.5 x 10 > n GR 4 c , M = m c + m x (8) 

Here, R c is the radius of the companion. Note that the strength of magnetic braking is still 
a matter of some controversy; while the evidence cited in the above reference argues for a 
reduction from the original value, it can also be argued that the presence of the well-known 
"period gap" in the period distribution of CVs requires a strength comparable to the original 
one. We have adopted here a recent prescription which is reasonably simple and adequate 
for our purposes: our final results do not depend significantly on the strength of this process. 

Consider finally collisional hardening. As in dicated e arlier , it is a stochastic process, for 
whose continuous limit we use the prescription of I Shull I (Il979l ). as has been done previously 
in the literature (see BG06 for a discussion). According to this prescription, the rate of 
increase of orbital binding energy E of a compact binary due to collisional hardening is 
given in this limit by: 

m 3 f 

= A c a-y, A c = 18G f — (9) 

m c m x 

Here, mj is the mass of the stars in the static background representing the cluster. We shall 
use M pc~ 3 and km sec -1 as the units of p and v c respectively. In the above units, the 
value of 7 for Galactic globular clusters typically lie between ~ 10 3 and ~ 10 6 (BG06). The 
relation between E and J is: 

j IE mx\ 

J 2E + 2\m, mx) 1 1 ' 

and the angular momentum loss rate is related to the shrinkage rate of the orbit a, or 
hardening, as: 

a J rh r rhv , s 

- = 2--2— -2— 11 

a J m c mx 

The m c and mx terms on the right-hand side of Eqn. ffTTl) are nonzero during mass transfer 
in the XB phase. In the PXB phase, m c = % = 0, so that a is simply related to J as (see 
BG06 and references therein): 

a J . 
- = 2- 12 
a J 
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Using Eqns. ffTOl) and ([9]), we have in this case, 




J \ 1 / E \ A c „ rrv f 

2 J 2 m c mx 

The total rate of loss of orbital angular momentum due to the above three processes is: 
3tot{o) = ( - ) = Jew (a) + 3mb{o) + jc(a) (14) 

TOT 



2.3.2. Hardening in X-ray binary (XB) phase 

As mass transfer starts upon Roche lobe contact, its effect on the angular momentum 
balance in the XB must be taken into account, in the manner described below. Note first 
that, for the radius of the Roche-lobe Rl of the companion, we can use either the 1971 
Paczyhski approximation: 

R L /a = 0.462 (^) ^ , (15) 



.M) 

which holds for < m c /mx < 0.8, or the 1983 Eggleton approximation: 

r, , 0.49 . , . 

RL,a = 0.6 + g 2 /3 1n(l + g-V3)> 9 = ( 16 ) 

which holds for the entire range of values of the mass ratio q. Both approximations have 
been widely used in the literature, and they give essentially identical results for the mass 
ratios of interest here. We have used the Paczyhski approximation here for simplicity of 
calculation. 

At the Roche- lobe contact point, Rl must be equal to the companion radius, the value 
of which is R c ~ O.6-R for a compani on of m r = .6M M (see above), according to the 



mass-radius relation for low mass stars (IGhosh 1 120071 ). For mx = 1.4M , this translates 



into an orbital radius of = 1.94i? at Roche lobe contact, using Eqn. fTl5l) . After this, 
the companion continues to remain in Roche- lobe contac t as the binary shrinks further, and 



continues to transfer mass (I van den Heuvel Hl99lL Il992l ). In other words, we have 



R c = 0.46a ( — £ ) , (a<a L ) (17) 



.M. 

throughout the XB phase. During this phase, the binary is already narrow enough that the 
collisional hardening rate is quite negligible compared to those due to gravitational radiation 
and magnetic braking. 
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Since no significant mass loss is expected from the XB in this phase, we have 



-m x . 



Combining Eqns. (TTTT) . (fTTl) and (jl~8j) with a mass-radius relation for the companion of the 
form 



R r oc mi 



we find: 



jtot{a)a (s - |) 



£ _i_ 5 

2 t 6 



M~m c 



(19) 
(20) 



Here, jtot{.a) = Jgw ( a ) +Jmb(o) is the effective total rate of loss of angular momentum, since 
the collisional-hardening contributions are negligible, as explained above. 

For the low-mass main sequence companions that we consider here, s ~ 1. However, 
when the mass of t he companion becomes less than about 0.03M Q , it becomes degenerate, 
so that s ~ —1/3 (IGhosh 1 120071 ). This results in a widening of the orbit (a > 0) from 
this point onwards, which we do not follow here, since our study is not aimed at such 
systems, as explained in Sec. 14.11 This change-over point is, of course, that corresponding 
to the well-known period mini mum of ~ 80 minutes in the orbital evolution of CVs and 
LMXBsf l van den Heuvel Ill992l ). Henceforth, we denote the value of a at the period minimum 
by a pm , and we terminate the distributions of a and n(a, t) in a at a minimum value of a pm 
in the figures shown in this work. Thus, in Fig. [21 we display the hardening rate a against a, 
beginning from a wide PXB phase, going into Roche lobe contact, and continuing through 
the mass-transfer XB phase upto the above period minimum. Note that a has a very weak 
dependence on a during the XB phase, which may have interesting consequences, as we shall 
see later. 



2.4. Binary formation processes 

Compact binaries with degenerate primaries and low-mass companions are formed in 
globular cluster (henceforth GC) cores primarily by means of two dynamical processes, 
namely, (i) tidal capture (tc) of a degenerate, compact star (white dwarf or neutron star) by 
an ordinary star, and (ii) an exchange encounter (exl) between such a compact star and a 
binary of two ordinary stars, wherein the compact star replaces one of the binary members. 
Accordingly, the total rate of formation of compact binaries per unit binary radius, R(a), 
consists of the above tc rate r tc (a) and exl rate r ex i(a): 



R(a) = r tc (a) +r exl (a) 



(21) 
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where a is the orbital radius of the compact binary so formed. We now consider the rates of 
formation by tidal capture and by exchange. 



2.4-1- Tidal capture 

In a close encounter between a compact star of mass mx and an ordinary star of mass 
m c with a distance of closest approach r p , tidal capture can occur if their relative speed v is 
less than an appropriate critical speed i>o(r p ), which we discuss belo w. The cross se ction for 



encounters within this distance r p is given by the well-known form (I Spitzer 1119871 ): 



, 7 2ixGMr v . 

°'=( m i + —jr- t ) ( 22 ) 



which gives the differential cross section for tidal capture around r p as: 

datc = | (2nr p + ^)dr p v < v (r p ) 
dr p \ v > v (r p ) 

The first terms in the right-hand sides of Eqs. ff22l) and fl23|) are the obvious geometrical 
cross sections and the second terms are due to gravitational focusing (also see below). It is 
clear that the latter terms dominate when r p is small, as is the case for the range of values 
of r p relevant to the problem we study here. We shall return later to the actual numerical 
values of r p of interest to us in this study. 

After being tidally formed, the binary is believed to circul arize very rapid ly to an orbital 



radius a = 2r p , assuming conservation of angular momentum (I Spitzer 1119871 ). Accordingly, 
the differential cross-section in terms of a is given by: 

d*tc = f + v<v (a) 

da \ v > v (a) 1 ' 



Here, vq{o) is the critical velocity in terms of a, obtained by setting r p = a/2 in Eq. ( 1251) 
below. 

In a sense, the whole cross-section as expressed above may be regarded as "geometrical", 
if we look upon pure considerations of Newtonian gravity as being geometrical. Details 
of the essential astrophysics enter only when we calculate the critical speed v (r p ), and 
an inversion of this relation (together with other plausible requirements; see below) then 
readily gives us the range of r p over which tidal capture is physically admissible. This is an 
interesting topic, with literature going back to the mid-1970s and earlier, and we summarize 
in this section those essential points which we need in this work. The basic physics of tidal 
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capture is of course that, during a close encounter, the degenerate compact star excites non- 
radial oscillation modes in the normal companion star through tidal forcing (in an encounter 
between two normal stars, each excites oscillations in the other): the energy required to 
excite these oscillations comes from the kinetic energy of relative motion of the two stars, so 
that if enough energy is extracted from this source by exciting these modes, the stars become 
bound after the encounter. This energy condition readily translates into one between vq and 
r p , giving an upper limit v on velocity for a specified r p as above, or, as expressed more 
commonly, an upper limit on the distance of closest approach r p for a specified velocity 
(actually, often a distribution of velocities, e.g., a Maxwellian, with a specified parameter in 
practical situations, as we shall see below). 

The above relation between Vq and r„ has been calculated in t he lite rature at various 



levels of detail. The pioneering estimates given in iFabian et.al. I (119751 ) or earlier works 
basically employ the impulse approximation for calculating the gain in the internal energy of 
the tidally-perturbed star, wherein the changes in the positions of the two stars dur i ng th e 



tidal interaction are neglected. A clear account of the procedure is given in I Spitzer I (119871 ). 
where the final result is evaluated for two normal stars of equal masses. Upon generalizing 
this procedure appropriately to the problem we study, where we have (a) unequal stellar 
masses mx and m c , and (b) the fact that only the normal star of mass m c undergoes tidally- 
induced oscillations, we obtain the following relation between v and r p : 

Mr P ) = fernxJ^j 2 r p l (25) 

Here, R m is the root-mean-square radius of the companion star, i.e., its radius of gyration 
which is gi ven in the polyr topic approximation as R 2 m j R\ ~ 0.114 in terms of the companion's 



radius R c (I Spitzer Ill987l ). 

To obtain the overall rate of tidal capture in the GC core of volume 47rrj?/3 per unit 
interval in a around a, we first consider this rate around a particular value v of the above 
relative velocity of encounter, i.e., r tc (a, v) = (4-7T /3)r^.kx p 2 \da tc / da)v , in terms of the above 
differential cross-section, remembering that the rate of encounter scales with the product of 
the densities kxp and p of compact stars and normal stars respectively. We then average 
this rate over the distribution of v, obtaining the form: 

4 

rtc(a) = -nrlk x p 2 (o- tc (a, v)v), (26) 

where the angular brackets indicate an average over the f-distribution. 

For the actual averaging, we adopt in this work a Maxwellian distribution f mx (v), as 
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has been widely done in the literature. A normalized Maxwellian is 

f mx (v) = Av 2 ex P (-Pv 2 ), /3=r^, A=-^=/3$, (27) 

2vi \/7i 



where v c is the velocity dispersion introduced earlier, for which we adopt the canonical value 
10 km s _1 in the numerical calculations (also see below). 

With the aid of Eqns. (jUj) , (l25|) and fl27|) . we perform the averaging and obtain: 

(a tc (a, v)v) = I geo + I grav , 
where, 

I geo = yfia[l- ex P (-(3v 2 (a))((3v 2 (a) + 1)] (28) 
I grav = 2^GMfi [1 - exp(-Pv 2 (a))] 

The terms I geo and I grav above arise due to what we described respectively as the geometrical 
term and the gravitational focusing term in the discussion below Eq. (|2"31) . Eqns. f[2T)|) and 
together give the total tidal capture rate as: 



'327T 3 

r lc u, ) i/ —kxTGM [1 - exp(-/^ 2 (a))] , (29) 

where T is the Verbunt parameter describing the total two-body encounter rate in the cluster 
core, as introduced earlier, and we have ignored I geo compared to Ig rav , which is an excellent 
approximation for the range of r p or a relevant here. 

We show in Fig. [3]r te given by Eq. (129]) as a function of a: this tidal capture cross-section 
is nearly constant for a < 5R Q , and decreases rapidly at larger a. At this point, we need to 
invoke additional physical arguments in order to estimate the range of values of a or r p over 
which tidal capture is actually possible, and use the above cross-section only over this range 
for our calculations. The lower bound to the above range comes from the requirement that 
the two stars must form a binary and not merge into each other, and the upper bound comes 
from the requirement introduced earlier that enough energy of relative motion between the 
two stars must be absorbed by the tidally-excited oscillation modes that the stars become 
bound. Consider the lower bound on r p first. Clearly, a minimum value of this bound must 
be the sum of the stellar radii, which in our case leads to the bound r p > R c rs 0.6i? Q . A 
more conservative bound comes from the requirement that the companion must underfill its 
Roche lobe after the binary has formed, i.e., R c < Rl, which, with the aid of Eq. (ITS]) and 
a = 2r p , yields r p > 1.6R C ~ R Q for the masses mx = 1.4M and m c = 0.6M Q we have 
here. The idea behind the latter requirement is apparently that if the companion overfills its 
Roche lobe at this point, the ensuing mass transfer is likely to lead to a merger. This seems 
reasonable at first, but detailed N-body simulations of recent years have suggested that this 



-16- 



require ment may, in fact, be too restrictive. In the simulations of 1 Portegies Zwart et.al. 
( Il997bl ). which included stellar evolutionary effects according to the scheme of these authors, 
systems which violated the latter requirement but satisfied the former one were allowed to 
evolve, with the result that details of the evolution determined which systems merged and 
which did not. In fact, these authors found a lower limit on a = 2r p of approximately 
a > Rq for tidal capture with an average companion mass very similar to ours, which is 
to be compared with the limits a > 1.2R Q from the first requirement above, and a > 2R & 
from the latter. In view of this, we have adopted the lower bound of a mm m 1.2R Q for our 
calculations here, as shown in Fig. [31 

Consider now the upper bound on r p . We have already given the relation between r„ 



and Vq by Eq. (|25p in the impulse approximation. Remembering that Vq = 1/(3 



for a Maxwellian, the above relation yields, for a canonical value v c = 10 km s as given 
above, an upper limit of r p < 10.2R C for a polytropic index n = 3 and one of r p < 14.1i? c 
for n = 1.5. Note that these bounds of r p /R c are larger than those given for t wo st ars of 
equal mass (roughly 8 for n = 3 and 11 for n = 1.5) in Table 6.2 of I Spitzer I (119871 ) by a 
factor of (mx/mc) 1 / 3 since r™^ / R c scales with the mass-ratio in this manner in the impulse 
approximation, as can be seen readily from Eq. (12"5]) . remembering that R c oc m c for the 
companions we consider here. That r™" 1 * /R c should increase with increasing (mx/m c ) is 
qualitatively quite obvious, since, other things being equal, a higher value of the mass ratio 
excites tidally-forced oscillations of larger amplitude. We return below to the question of 
the exact scaling with this mass ratio. 

As has been realized long ago, the impulse approximation is of limited validity, work- 
ing best when the frequency of perturbation (i.e., tidal forcing) is not very different from 



those of the stellar oscillation modes that are excited by this perturbation (IFabian et.al. 



19751 ; I Spitzer 1119871 ). Since this is not the case for the values of r™ aa: /i? c estimated above, 
we need more accurate results, which come from detailed computations of the total en- 
er gy dissipated by the above excited modes. Such numerical computations were pioneered 
by I Press &: Teukolsky I (119771 ). and detailed results were established for vari ous situations by 
severa l groups of authors in the mid-1980s, including I Lee &: Ostriker I (119861 ) and l McMillan et.al. 
( 119871 ). which have been extensively used since. These results have shown that the exact 
upper bounds on r p are considerably smaller than those given by the impulse approximation, 
as may have been expected, since the forcing frequency falls far below those of the oscillation 
modes at such large separations as are given by this approximation, and the effi ciency of 



excitin g these modes drops rapidly. Some exact results are given in Table 6.2 of I Spitzer 



( 119871 ) from the above references, but only for the equal-mass case, where the above upper 
/ R c is 2.4 for n = 3 and 3.4 for n = 1.5. 



bound r™ aa 
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For our purposes here, we need to obtain the above upper bounds for our mass ratio 
mx/m c = 1.4/0.6 ~ 2.3, which we do by doing a power-law fit of the form r™ ax oc (mx/m c ) a 



to the results given 



or va rious values of the degenerate/normal star mass-ratios in Table 3 



of I Lee fc Ostriker I (119861) . This yields a ~ 0.62 (note that the quantity listed in Table 3 of 
Lee fc Ostriker ( 19861 ) is the impact parameter R defined by these authors; r™ ax scales as 



Rq, as shown in their paper). The interesting point about this scaling is that it is stronger 
than that given above by the impulse approximation, which corresponds to a = 1/3. Clearly, 
then, the impulse approximation fails to extract the entire scaling with mx/m c . The reason 
for this appears to be related to nonlinear effects in exciting and dissipating tidally-induced 
oscillations, but needs to be investigated furtheiB With the above value of a, the upper 
bound r™ ax ' /R c for our mass-ratio here is 4.1 for n = 3 and 5.7 for n — 1.5. As the latter 
value of the polytropic index is believed to give a better representation of a low-mass main- 
sequence companion of the kind we are considering here, we adopt r™ ax /R c w 5.7 here. 
With a = 2r p and the value of R c given earlier, this translates into an upper bound on a as 
a max ^ g,3^ Q; w hich we can adopt for these calculations. 

Thus we find a range of values 1.2i2 < a < 6.8R Q over which tidal capture is expected 
to be effective in the problem we study here. Consider now how the tidal-capture cross- 
section is expected to fall off at the bounds of this range. At the upper bound, the cut-off 
is not sharp, of course, as there is a distribution of velocities. In other words, the upper 
bound a max as given above corresponds to a suitable average (actually, root-mean-square 
in this case) velocity, so that at any a > a max , there will be some stars in the distribution 
whose velocities are sufficiently below this average that tidal capture will be possible for 
them. Of course, their number will decrease as a increases, producing a "tail" in the tidal 
capture cross-section whose shape is determined by that of the velocity distribution. We 
have used a Maxwellian distribution here, which gives the tail seen in Fig. [3j which falls off 
rapidly beyond a max = 6.8i? Q . We shall use this fall-off profile in our calculations: other 
profiles will not make a large difference. At the lower bound, in view of the discussion 
given earlier, we expect the cross-section to actually fall off gradually from about a = 2R Q 
to a = a mm = 1.2R & , rather than being cut off sharply at a mm , but we shall ignore this 
complication here. 

We close this discussion of tidal capture with some observations on the many inves- 
tigations, conclusions, and points of view that the subject has now seen for more than 



2 Note that this discrepancy is even stronger for the case where both stars are normal, main-sequence ones, 
since a ~ 1.6 in that case, as can be shown readily from Table 2 in the above Lee-Ostriker reference. An 
obvious line of reasoning for this would be that larger nonlinear effects may be expected when two normal 
stars force tidal oscillations in each other, but we shall not speculate on this any further here. 
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three decades . From the pioneering suggestion and an essentially dimensional estimate of 



Fabian et.al. I (119751 ). detailed calculations of the 1980s and '90s have reached interesting, 



and sometimes contradictory, conclusions. For example, concerns that energy dissipation by 
tidally- induced modes may lead to a large distention of the companion and so to a merger 
have been confronted with results from detailed computations of the nonlinear damping 
of the primary modes by coupling to other, high-degree modes, which suggested that the 
damping took place far more rapidly than thought before, and the energy dissipated was too 
small to have a significant effect on the companion's structure. We here have a adopted a 
somewhat moderate view that tidal capture is plausible, but efficient over only a restricted 
range of r p or a. This view is supported by (a) recent observational demonstration that the 
number of X-ray sources in Galact ic globular clusters s cale with their Verbunt parameter T, 
i.e., the two-body enc ounter rate (jPoolev et al. 1120031) . as described earlier, and (b) recent 



N-body simulations of I Portegies Zwart et.al. I (ll997bl ) showing tidal capture over a consid- 



erable range of a, admittedly under the algorithms adopted by these authors. Consider, 
finally, our suggested range of radii for efficient tidal capture, a max /a mm « 5.7 
above, in the context of other suggested ranges. Va l ues in the ra nge a max /d 



as given 

imil ~ 2 3 



have be en thought plausible by I Podsiadlowski et.al. I (120021 ) , while I Portegies Zwart et.al. 
(1l997bl ) have demonstrated tidal capture over a range a max /a mm ~ 10. We here advocate a 
range a max j a mm ~ 4 — 6 (depending on n), which is between the two, and still quite modest. 



2.4-2. Formation by exchange 



Exchange encounters be tween binaries and single s t ars wi th arbitrary mass ratios has 
been extensively studied by I Heggie. Hut fc McMillan I (119961 ). They performed detailed 
numerical scattering experiments, using the automatic scattering tools of the STARLAB 
package. From the resulting exchange cross sections, they obtained a semi-analytic fit of the 
form: 

nGM tot R_ 



2v 2 



-a{m l ,m 2l m 3 ) 



(30) 



Here, R is the orbital radius of the initial binary, mi is the mass of the escaping star, m 2 
is the companion mass, is the mass of the incoming star, and M to t = mi + mi + m 3 . 
(j{ra\, 7T7.2, rn.3) is the dimensi qnless cross section which is a func tion of these masses only and 
which is given by Eq. (17) of I Heggie. Hut Sz McMillan I (ll9961 ). We use Eqn. f[3"Uj) to obtain 



the cross sections a ex i(a) for the exchange process 'exl' described above, but one essential 
point needs to be clarified first. 



The radius a of the compact binary formed by exchange is not the same as the radius 
a' of the original binary undergoing exchange. Therefore, a relation between a' and a is 
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required, since in Eqn. fl30l) R represents the radius a' of the initial binary, not the radius a 
of the compact binary formed by exchange. According to the binary-hardening rule of Heggie 



( IHeggie 1119751 ). the final compact binary must, on an average, be harder, i.e., have a larger 
binding energy. We performed illustrative scattering experiments with circular binaries and 
incoming stars with mass ratios of interest to us in this study, using the scattering tools 
of STARLAB. The resulting distribution of the change in orbital radius Aa/a is shown in 
Fig. HI and is seen to be highly asymmetric. 

The long tail towards Aa > implies that the binary radius increases in many scatter- 
ings. This does not of course contradict the above Heggie rule, since the increase of mass 
due to exchange (the mass of the incoming compact star, 1.4M , is a factor « 2.3 times the 
mass of the outgoing low-mass star, O.6M ) increases the binding energy by itself by the 
above factor. From these experiments, we see that the peak of the distribution corresponds 
to a shrinkage of the binary by about 25 per cent. On the other hand, the average change 
in binary radius, calculated from the above distribution, is much closer to zero due to the 
above long tail of the distribution on the Aa > side, so that we can take a m a' for our 
purposes here without much error. 

The total Maxwellian-averaged rate of formation of compact binary by this type of 
exchange (exl) in the GC core is then: 



4 8n 3 
r«Bi(a) = -nrlkxp 2 fb(a){o- ex i(a)v) = d —kxfb(a)TGM tot aW(m c ,m x ) (31) 

Here, /&(a) is the distribution function of the orbital radii of the primordial stellar binaries 
in the cluster core. For primordial binaries, we can take the widely- used distribution /&(a) oc 



1/a (i.e., a uniform distribution in In a) (I Kraicheva et.al. 1119781 ). with a lower bound at 



a fa 1.2R & , corresponding to the smallest possible radius for a binary of two 0.6M Q main- 
sequence stars. The exl rate is shown in Fig. [3j 



2.5. Binary destruction processes 

A compact binary can be destroyed by two major processes. First, an encounter with a 



star which has a relative speed higher than an appropriate critical speed (I Hut &: Bahcall 



19831 ) can lead to its dissociation (dss). Second, in an exchange encounter (ex2) of this binary 
with a compact star, the latter can replace the low-mass companion in the binary, forming a 
double compact-star binary consisting of two neutron stars, two white dwarfs, or a neutron 
star and a white dwarf, all with masses mx ~ 1.4M . This, in effect, destroys the binary 
as an X-ray source (as accretion is not possible in such a system) , and so takes it out of our 
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reckoning in this study. This is so because such a system is not an X-ray source, and it is 
essentially impossible for one of the compact stars in such a system to be exchanged with 
an ordinary star in a subsequent exchange encounter, since rrif = O.QM & is much lighter 
than m x = 1.4M . The total destruction rate D(a) per binary is thus the sum of the above 
dissociation and exchange rates: 

D(a) = r ex2 {a) + r dss (a) (32) 
We now discuss the rates of these two processes. 



2.5.1. Dissociation 

To estima te the dissociation rate of compact binaries, we use the results of scattering 



experiments of I Hut &: Bahcall I (I1983I ). The Maxwellian-averaged dissociation rate (dss) per 



compact binary is then given by 

rdss{a) = k x p{a dss (a)v) (33) 



From I Hut fc Bahcall I (119831 ). we adopt 



(a dss (a)v) = ^J-v c a 2 exp . (34 ) 



32tt 1 6 2 / 3v 2 
— ^-v c a exp (--- ? 

a relation which was obtained by these authors by fitting the results of their scattering 
experiments with analytical models. Here, v cr u is the threshold relative velocity for ionization 
(see Sec. I2.5p . given by: 

2 _ Gm x (2m c + m x ) 1 

V crit ~ i • {^d) 



As these authors pointed out, Eqn. ( |34j) is an asymptotic form, which works well only for 
significantly hard binaries, i.e., those with v c « v cr n. This condition is of course satisfied 
for the compact binaries that we are interested in here. 

We show in Fig. [3] the above dissociation rate, whose essential variation with a is seen 
by combining Eqs. ( 134"1) and ( 1351) . which yields the form rd ss (a) oc a 2 exp(— a c /a), where a c 
is a constant. Thus, the dissociation rate is quite negligible for aCa C) reflecting the fact 
that it is essentially impossible to dissociate very hard binaries. As a increases, the rate rises 
extremely sharply at first (the initial rise is determined by the exponential), and eventually 
scales as a 2 for a ^> a c . 



2.5.2. Destruction by exchange 



By arguments similar to those given in Sec. 12.4.21 we arrive at a Maxwellian-averaged 
rate of this type of exchange (ex2) per compact binary which is: 

r e x2(a) = k x p(o- ex2 (a)v) = d —k x ^GM tot aW(m c ,m x ), (36) 

and which is also shown in Fig. [3j This rate scales with a simply as r eX 2(a) oc a. Note the 
different magnifications used for different curves in Fig. [3] in order to make all of them clearly 
visible. Of the two destruction processes, r ex2 dominates completely at all orbital radii of 
interest in our study (reflecting the fact that dynamically-formed binaries in GC cores are 
so hard that they cannot be dissociated or "ionized" by further encounters in that GC core), 
but the fast-rising r<i ss eventually overtakes it at a ~ lOOOi? , corresponding to very soft 
binaries. 



2.6. The numerical method 

Equation (|5]) for the evolution of compact binary populations is a partial differential 
equation (PDE) of hyperbolic type, with similarities to wave equations, as point ed out 



earlier. We solved this equation using a Lax-Wendorff scheme (I Press et.al. Ill992l ). This 



involves dividing the range of a and t in a discrete mesh (aj,tjv) of constant space intervals 
(Aa) and time intervals {At). The PDE is then discretised into a set of linear difference 
equations over this mesh, which is solved numerically. 

We denote by nf the value of n at the iVth time step and the jth point in a. Discreti- 
sation of Eqn. ([5]) according to the Lax-Wendorff scheme is a two-step process: 



Half step : 

AT+l/2 
l j+l/2 - 2 
/0j+i/2)Af , N 



a . , \ » = \ (nf +l + nf) + R(a j+1/2 ) - D(a j+1/2 ) 



2Aa 



in 



'j'+i 



Full step : 



nf +1 = nf + (R{ aj ) - D{ aj )nf) At 



/(°j)At (N+l/2 _ N+l/2 
Aa \ n j+l/2 n j-l/2 



n N +n N 



At 
2 



(37) 



For a chosen mesh-interval Aa, Eqn. fl37|) will be numerically stable o nly if the time-step 
At is chosen to be small enough that it obeys the Courant condition (I Press et.al. Ill992l ) 
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throughout the mesh: 



At = 7] 



Aa 

fmax 



rj < 1 



(38) 



where f ma x is the maximum value of f(a) within the a-range of the mesh. 



We chose Lax-Wendorff scheme among the various existing schemes for solving hyper- 
bolic PDEs primarily because it appe ars to be the only e xplicit method that does not have 



any significant numerical dissipation (1 Press et.al. Ill992l . and references therein) and is at 



the same time numerically stable, provided that the time step is chosen according to the 
Courant condition. This point is important, since numerical dissipation can significantly af- 
fect the computed evolution of n(a, t) and the X-ray binary population, as we observed while 
trying other methods, e.g., the so-called staggered -leapfrog method. Other instabilities, e.g., 



the mesh-drifting instability (1 Press et.al. 
method we chose. 



19921 ). also appeared to be insignificant in the 



3. Results 

3.1. Evolution of compact-binary distribution 

A typical result from our computed evolution of the compact-binary distribution func- 
tion n(a, t) is shown in Fig. [SJ wherein the surface n(a, t) is explicitly displayed in three 
dimensions. The GC parameters chosen for this run were p = 6.4 x 10 4 M & pc~ 3 , r c = 0.5 pc 
and v c = 11.6 km sec -1 , similar to those of the well-known Galactic cluster 47 Tuc. The 
distribution function is seen to evolve as a smooth surface, with the compact binary popula- 
tion growing predominantly at shorter radii (a < 10i? Q , say). We start with a small number 
of binaries at t — following various distributions, and find that the distribution at large 
times ~ Gyr is quite independent of these initial conditions, being determined entirely by the 
dynamical processes of formation and destruction, and by the various hardening processes 
detailed earlier. Note that, since the point of Roche lobe contact corresponds to a « 2R Q 
in our study, as explained earlier, that part of the distribution which is shortward of this 
radius corresponds to XBs, while that part longward of it corresponds to PXBs. 

To further clarify the nature of this evolution, slices through the above surface at various 
points along time axis and a-axis are shown in Figs. [6] and [71 in the former figure the abscissa 
being also marked in terms of the orbital period P, readily calculable in terms of a and the 
stellar masses with the aid of Kepler's third law, assuming conservative mass transfer during 
the XB phase. Figure [6] shows that n(a) increases with time, roughly preserving its profile 
for t > 1.5 Gyr or so. This profile consists of a roughly uniform distribution in for short 
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orbital radii, a < QRq, say, corresponding to P < l d roughly, and a sharp fall-off at larger 
radii and orbital periods. Figure [7J shows that n(a) at a given a increases with time and 
approaches saturation on a timescale 6 — 12 Gyr or so, this timescale being longer at at 
smaller values of a. 

Figures [6] and [7] suggest that a regime of roughly self-similar evolution may be occurring 
in our model binary population at times beyond 1 Gyr or so, in the following way. An 
asymptotic profile of n(a) is established on the timescale of a 1 Gyr or so, which thereafter 
evolves roughly self-similar ly towards a saturation strength on a timescale ~ 6 — 12 Gyr 
or so. We shall discuss the origins of such behavior in detail elsewhere, since, as explained 
in Sec. 14.11 our model of orbital evolution requires additional ingredients before it can be 
compared with observations of X-ray binaries. However, the following qualitative remarks 
are appropriate here. 

First, the origins of the establishment of the above self-similar profile in a Gyr or so 
(independent of the initial distribution we start from) clearly lie in the two terms that 
describe binary formation and hardening on the right-hand side of Eq. (151), namely, R(a) 
and |p/(a) respectively. The latter term can be written qualitatively in the form n/rh, 
where is the overall hardening timescale, which is well-known to be of the order of a 
Gyr or so (see BG06 and references therein). This timescale, which is also that on which a 
given binary passes from the large-a end of the distribution shown in Fig. [5] to the small-a 
end, is obviously the timescale that establishes the above profile. The shape of this profile, 
as detailed above, seems related to those of the tidal-capture rate (see Fig. [3]) and the 
hardening rate (see Fig. [2]). In particular, note that the former rate is roughly constant over 
a min < o < 5-R Q , and the latter roughly so for a pm < a < 2R Q . 

Second, the subsequent, roughly self-similar evolution of the above profile occurs on a 
(longer) timescale t s whose origins clearly lie in the binary destruction term on the right- 
hand side of Eq. (Q, namely, nD(a), since this term can be cast in the qualitative form n/r s , 
where r s is the saturation time ~ 6 — 12 Gyr. Whereas the earlier term n/rh describes the 
passage or "current" of binaries through the distribution, as described earlier, the term n/r s 
becomes important as n increases, preventing n from becoming arbitrarily large by enforcing 
saturation at the point where the rates of formation and destruction balance. As D(a) scales 
with a, as shown above, and t s = 1/D(a), we expect saturation to occur at earlier times at 
larger radii, as seen in Fig. [7J 



-24- 



3.2. Number of X-ray binaries in globular clusters 

The total number of X-ray binaries Nxb i n a GC at any time can be computed directly 
from our approach by integrating n(a, t) over the range of a relevant for XBs, viz., a pm < a < 
cil, where a pm is the value of a corresponding to the period minimum P « 80 minutes, and 
a L is the value of a at the first Roche lobe contact and onset of mass transfer, as explained 
earlier. We have: 

N XB {t) = / n(a,t)da (39) 

Taking an evolutionary time ~ 8 Gyr as representative, we can therefore determine Nxb 
at this point in time, and study its dependence on the Verbunt parameters T and 7 that 
describe the essential dynamical properties of globular clusters in this context, as explained 
earlier. By doing so, we can attempt to make qualitative contact with the systematics of 
those r ecent observations o f X-ray binaries in globular clusters which we have described 



earlier (jPooley et al. 1 120031 ) . To this end, we computed values of Nxb over a rectangular 
grid spanning over 7=1 — 10 6 and T = 10 3 — 10 8 . (Of course, not all the points on the 
grid would be directly relevant for comparison with observation, since the observed globular 
clusters lie only along a diagonal patch on this grid, as shown in Fig. [IJ However, in this 
introductory study, we wished to establish the theoretically expected trends of variation with 
T and 7, and so performed computations of Nxb over the entire rectangular grid) 

For a specified grid point, i.e., a pair of values of the Verbunt parameters, we obtained 
values of p, r c and v c with the aid of Eqs. (JTJ, (j2J) and the virialization condition: 

v c oc p 1/2 r c (40) 

which were used for the computation at this grid point. We chose this prescription for the 
sake of definiteness, because values of v c are not known, in general, at a computational grid 
point, without which a pair of Verbunt parameters cannot specify all three variables p, r c 
and v c . This also introduced a certain uniformity of treatment of all grid points, which, we 
thought, would clarify the theoretically expected trends. On the other hand, this did lead 
to a feature at high values of F and low values of 7, i.e., in that part of the grid which is 
completely devoid of observed globular clusters at this time (and which, in fact, may actually 
contain no clusters, because such combinations of T and 7 may not be possible in nature), 
which appears unphysical, as we discuss below. Observationally, we know, of course, that 
some clusters appear fairly virialized and some do not, but any spread in v c applied over 
the grid points would have been arbitrary, and would have led to a scatter, masking the 
systematic theoretical trends without purpose. Finally, throughout these computations, we 
used representative values for (a) the primordial binary fraction kb, namely, 10 percent, and 
(b) compact star fraction kx, namely, 5 percent. 
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Figure [8] shows the computed surface Nxb(j,T). There appears to be a "fold" in this 
surface, in a direction roughly parallel to the F axis, located around 7 = 3 x 10 3 . From this 
fold, if we go towards higher values of 7, then, for any given value of T, N X b decreases with 
increasing 7. This is a signature of the compact-binary destruction processes detailed in the 
previous section, whose strengths increase with increasing 7. Thus, the above value of 7 
corresponding to the fold seems to be a good estimate of the threshold above which these 
destruction processes dominate. At constant 7, the variation with V is quite straightforward: 
Nxb simply increases monotonically with increasing T, reflecting the fact that the formation 
rates of compact binaries, as described in the previous section, increase with increasing T. 

To further clarify these trends, and to facilitate comparison with those obtained from 
the "toy" model of BG06, we display in Fig. V /Nxb vs. 7, as was done in that reference. 
The motivation is as follows. It was shown in BG06 that the toy model of these authors leads 
to the scaling that V /Nxb was a function of 7 alone, which was a monotonically increasing 
function of 7, for which the toy model gave a very simple, analytic form. Our purpose in 
Fig. [9] is to see how much of this scaling survives the scrutiny of a more detailed model, 
such as presented here. As the figure shows, this scaling does carry over approximately, 
although some details are different. T/N X b is still almost a function of 7 alone (except 
at the very highest values of T), showing that this scaling N X b oc Tg(j) of the toy model 
carries over approximately to more detailed ones, thereby giving an indication of the basic 
ways in which dynamical binary formation and destruction processes work. The above 
"universal" function (7(7) of 7 is, except for a feature at low values of 7 which we discuss 
below, still a monotonically increasing one, reflecting the increasing strength of dynamical 
binary-destruction processes with increasing 7. However, the shape of the function is different 
in detail now, as may have been expected. 

We now discuss the low-7 feature referred to above: at the lowest values of 7, T/N X b 
seems to rise again, reflecting an apparent drop in N X b- This is difficult to understand, since 
binary-destruction effects are negligible at these values of 7. Actually, this is an artifact of 
the way in which we fixed the essential cluster parameters p, r c and v c from specified values 
of the Verbunt parameters for the computational grid (as explained above), which can be 
seen as follows. With the assumption of virialization, as done for this purpose, the velocity 
dispersion v c can be expressed in terms of the Verbunt parameters in a manner analogous to 
that used in Eq. ([3]), the result being v c oc r 2 / 5 7 _1//5 . This relates v c to 7, so that the latter 
influences the Maxwellian-averaging process involved in the calculation of the tidal capture 
cross-section described in Sec. 12.4.11 since the parameter (3 = 3/(2t> 2 ) of the Maxwellian 
then scales as (3 oc r~ 4 / 5 7 2 / 5 . At small values of 7, (3 becomes small, which reduces the 
tidal-capture rate, as Eq. fl29|) readily shows. This is completely unphysical, of course, since 
7 has nothing to do physically with the tidal capture rate. Rather, it is an artifact produced 
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by the way we (artificially) related v c to 7 for computational convenience. Accordingly, we 
ignore this low-7 feature in all further considerations. 



4. Comparison with Observation 
4.1. Applicability of our study 



Before attempting to compare our results with observations, we review in brief some 
essential ingredients of our model study at this stage, so as to clarify which of our results 
can be so compared, and which need inclusion of further components before this can be 
meaningfully done. A major ingredient that is incomplete at this stage is our description 
of the orbital evolution of the binary, since it neglects nuclear evolution of the low-mass 
companion star altogether. While this may not be very unreasonable for CV systems or 
for short-period LMXBs with orbital periods between ~ 10 hours and the above period 
minimum of ~ 80 minutes, it is completely inadequate for other LMXBs, where the stellar 
evolution of the companion pla ys a crucial role, whic h has been studied by many au thors. 
In particular, recent studies by I Podsiadlowski et.al. I (120021 ) and I Pfahl et.al. I (120031 ) have 
demonstrated the large range of possibilities covered by such evolution with realistic stellar 
evolutionary codes, performing a Monte Carlo binary population synthesis study in the latter 
reference with the aid of the library of evolutionary sequences described in the former. We 
plan to include stellar evolutionary effects in a subsequent work of the series and are assessing 
various methods of doing so. One possibil ity is to start with a semi-analy tic scheme along 
the lines of the SeBa model as described in I Portegies Zwart et.al. I (Il997bl ). and to continue 
with a semi-analytic approximation to a more elaborate library of evolutionary sequences, 
such as described above. 



Since most of the XBs in the Galactic GC data of iPooley et al. I (120031 ) are CVs, our 
scheme should be able to describe the overall properties of these XB populations reasonably 
well. Even so, we shall make no attempt here to compare our results on orbital period distri- 
bution with the observed CV distribution, since the CVs in the latter distribution are almost 
exclusively from outside globular clusters, where dynamical formation is not relevant. We 
have here recorded the orbital-period distribution that comes from our computations (at this 
stage) only as a natural intermediate step. It can perhaps be compared with observation 
when the orbital-period distribution of CVs in GCs becomes observationally established. 
For LMXBs, where the observed orbital-period distribution at this time also consists over- 
whelmingly of those outside GCs, there is of course no question of comparison at this stage, 
for the reasons given above. Thus, our main aim here is to put in the observational context 
our results on the numerical properties of XB populations in GCs in relation to the GC 
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parameters. 



4.2. Ultracompact X-ray binaries 



In recent years, a subset of LMXBs in GCs, in the Milky Way and possibly also in ellip- 
tical galaxies, have received much attention because of (a) their high, persistent brightness 
(L x ~ 10 36 — 10 39 erg s _1 ), which would make them dominate the h igh end of the luminosity 



functions of X-ray binaries in ellipticals (I Bildsten fc Deloye 1120041 ) and (b) their very close 
orbits with P < 1 hr or so, sometimes as short as P ~ 10 minutes, the classic example being 
the 11 min binary 4U 1820-30 the Galactic cluster NGC 6624. These are the ultracompact 
X-ray binaries (henceforth UCXBs), which are thought to consist of neutron stars in ultra- 
compact orbits with very low-mass degenerate dwarf companions (m c ~ (0.06 — O.2)M ) as 
mass donors. The evolutionary origin of UCXBs is of much current interest, and proposals 
for such origins include (a) direct collisions between red giants and neutron stars in GC 
cores, as a consequen ce of which the red-giant envelope can either be promptly disrupted 
(1 Ivanova et.al. 1120051 ) or be expelled more slowly in a common-envelope phase, and (b) usual 
LM XB evolution with the initia l orbital period below the "bifurcation period" of about 18 
hrs (I Podsiadlowski et.al. II2002I ). A natural point that arises, therefore, is about the role 



of UCXBs in our study, and the general importance of the above channels of formation in 
relation to the ones we have described above, which we now consider in brief. 

The key feature of UCXBs from the point of view of our study is that the number of 
UCXBs Njjc is a tiny fraction of the total number of XBs in a GC, and so of little importance 
as far as Nxb is concerned. This is a general, robust feature, which follows from the basic 
point that the UCXBs are extremely short-lived because of their extreme brightness, so that 
Nuc is small at any given epoch despite their considerable birthrate. To see this in more 
detail, consider the UCXB birthrate of about one every 2 x 10 6 year per 1O 7 M of the mass of 
a GC, as given by (I Bildsten &: Deloye 112004 ). which, together with their estimated lifetimes 
of (3 — 10) x 10 6 years, yields an estimate of Nuc ~ 1 — 5 in a 1O 7 M GC at any given time. 
Actually, the observed GCs in our galaxy have lower masses, in the range ~ (10 5 — 1O 6 )M 
J Ivanova et.al. hood ). Thus a Galactic GC of 1O 6 M like 47 Tuc will have N uc ~ 0.1-0.5, 
remembering that the birthrate scales down appropriately with the GC mass, but the lifetime 
re mains the same. Thi s is to be compared with the observed number of XBs in 47 Tuc of 
45 JPoolev et al. Il2003h . which yields a fraction Nuc/N XB ~2x 10~ 3 - 1.1 x 10~ 2 . We can 



doubl e- check this and put it on a systematic basis with the aid of Table 1 of I Ivanova et.al. 



(120051 ). wherein these authors have listed the minimum expected number of UCXBs in a 
number of Galactic GCs, by combining this with the total number of observed XBs obtained 
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from iPooley et al. I (120031 ) and other sources. For 47 Tuc, with 0.23 UCXBs and 45 XBs, 
the ratio is Nye /Nxb ~5x 10~ 3 , very similar to above, and those for other sources are also 
similar. For example, Terzan 5 has a ratio ~ 2 x 10~ 3 , and NGC 6652 has a ratio ~ 8 x 10~ 4 . 

We see from the above that UCXBs constitute such a tiny fraction of the total XB 
populations of Galactic GCs in terms of numbers that their effect is negligible for this work. 
However, in a study of the X-ray luminosity functions of GCs, their effect is expected to be 
crucial: if a GC contains even one UCXB, its luminosity may dominate over the combined 
output of all other XBs. It is the extension of this idea which has been used in recent years 
to argue th at the luminosity function of XBs in ellipticals may be dominated by UCXBs in 
their GCs J Bildsten fc Delove lliooi l. 



4.3. X-ray source numbers in globular clusters 

The filled squares in Fig. [S] represent globular clusters with significant numbers of X-ray 
binaries in them. These points generally lie near the surface in this three-dimensional plot, 
mostly in the vicinity of the fold described above. This is more clearly seen in the two- 
dimensional plot of Fig. [9X where the bulk of the observational points are indeed seen to be 
near the upward "knee" of the computed curves. To facilitate comparison with observations 
further, we show in Fig. [TU] contours of constant Nxb in the T— 7 (Verbunt parameters) plane. 
Overplotted on these are the above observed clusters (filled sqaures), where the number in 
the parentheses nex t to each indicates the total number of X-ray binaries observed in it 



(IPooley et al. 1120031 ). The contours are seen to be qualitatively rather similar in shape to 
the curves in Fig. The trend in the observed Nxb values generally follows the contours, 
with one exception. This is most encouraging (also see BG06) for the construction of more 
detailed models, and indeed rather remarkable in view of the fact that no particular attempt 
has been made to fit the data at this stage. 



5. Discussion 

In this paper, we have explored the results of a Boltzmann study of the evolution 
of compact-binary populations in globular clusters in the continuous limit, and made pre- 
liminary contacts with observations of X-ray binaries in Galactic globular clusters. Our 
Boltzmann approach has built into it the rates of the essential dynamical processes that 
occur due to star-star and star-binary encounters in dense clusters, viz., collisional harden- 
ing, binary formation by tidal capture and exchange, and binary destruction by dissociation 



-29- 



and other mechanisms, as obtained by previous numerical studies of large numbers of such 
individual encounters. We stress that our Boltzmann scheme is not a Fokker-Planck one, 
wherein the cumulative effects of a large number of small changes in distant encounters is 
described as a slow diffusion in phase space. We can and do handle both small and large 
changes in the framework of the original Boltzmann visualization of motion through phase 
space (at a computational cost which is quite trivial compared to that required for N-body 
simulations). Indeed, the continuous limit of collisional hardening used in this paper may be 
looked upon as an example of a slow diffusion in a-space, while some of the formation and 
destruction processes are examples of faster and more radical changes. Of course, all these 
processes are episodic in nature, and we are studying their continuous, probabilistic limit in 
this introductory paper. As already pointed out, Paper II will describe an explicit treatment 
of the stochasticity of these processes within the framework of stochastic PDEs, which the 
Boltzmann equation becomes under such circumstances. 



5.1. Conclusions 

We find the indications from this preliminary study to be sufficiently encouraging to 
attempt several steps of improvement, most of which we have already indicated in the pre- 
vious sections. To recapitulate briefly, we need to provide an appropriate description of the 
stochastic processes, which we do in Paper II. We need to introduce a mass function for 
the background stars in the globular cluster core, and handle non-circular orbits formed in 
the encounter processes. We need to assess the possible importance of binary-binary interac- 
tions in this problem, which we have neglected altogether so far. We need to include essential 
aspects of stellar evolution of the companion in our orbital-evolution scheme, particularly 
for LMXBs. In a more ambitious vein, we need to consider the evolution of the stellar 
background representing the cluster core, which we do in Paper III. As the core collapses, 
the collapse stalls due to binary heating, and possible gravothermal oscillations occur, the 
core parameters p and r c evolve appropriately, and so do the Verbunt parameters T and 7. 
Binary-population evolution with such evolving GC parameters is an interesting problem 
in itself, even if we do not explicitly consider the back reaction of binary evolution on the 
evolution of its background. 

The scaling of Nxb with the two Verbunt parameters we already found here seems to be 
among the basic building blocks of our understanding of how globular clusters cook up their 
gross overabundance of X-ray binaries through an interplay between dynamical formation 
and destruction. It remains to be seen if there are other such building blocks which have not 
been investigated so far. 
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Fig. 1. — Positions of Galactic globular clusters (open circles with dots) on the T — 7 
(Verbunt parameters) plane. Overplotted are positions of those clusters with significant 
numbers of X-ray sources detected in them (filled squares), color-coded according to the 
number of X-ray sources in e ach, the color code being displayed on the right. Data from 



Harris I (11996. revised in 19991 ). 
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Fig. 2. — Hardening rate a of a compact binary as a function of the orbital radius a, in 
a globular cluster with a Verbunt parameter of 7 = 10 3 . Collisional hardening dominates 
roughly at a > 14i? , and gravitational radiation plus magnetic braking roughly in the 
range 2R Q < a < 1AR Q . These two regions, shown as dashed lines, are populated by pre- 
X-ray binaries (PXBs), which are detached. At a m 2R Q , Roche lobe contact occurs and 
mass transfer begins, so that the region shortward of this radius, shown as the solid line, is 
populated by X-ray binaries (XBs). This region is shown upto the orbital radius a pm which 
corresponds to the period minimum of ~ 80 min (see text). Along abscissa, both orbital 
radius a and orbital period P scales are shown for convenience. 
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Fig. 3. — Tidal capture (tc) rate, the exchange rates 'exl' and 'ex2', and the dissociation 
(dss) rate, as described in text. Note that, compared to the tc rate, the exl rate has been 
magnified by a factor of 50, the ex2 rate rate by a of factor 60, and the dss rate by a factor 
of 10 9 , so that all rates are clearly visible. Along abscissa, both orbital radius a and orbital 
period P scales are shown for convenience. Curves are terminated at a radius a m i n = 1-2R& 
(see text). 




Fig. 4. — Distribution of the fractional change in binary radius Aa/a for ~ 30000 scattering 
experiments with v/v cr u = 0.5 (see text) and random impact parameters. The distribution 
is highly asymmetric, with a peak at Aa/a ~ —0.25, and a long tail in the Aa > direction. 



-37- 




Fig. 5. — Three-dimensional surface n(a,t) describing the model evolution of population- 
distribution function of compact binaries for GC parameters p = 6.4 x 10 4 M pc~ 3 , r c = 
0.5 pc, v c = 11.6 km sec -1 (roughly corresponding to 47 Tuc). The lines on the surface 
represent only samples from the set of computed points, the computation having been done 
over a much finer grid. 
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Fig. 6. — Time slices, i.e., n(a) at specified times t, for the evolution n(a,t) shown in Fig. [5] 
Along abscissa, both orbital radius a and orbital period P scales are shown for convenience. 
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Fig. 7. — Radial slices, i.e., n(t) at specified orbital radii a, for the evolution n(a,t) shown 
in Fig. O 
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Fig. 8. — Computed Nxb(T, 7) surface. Overplotted are the positions of the globular clusters 
with significant numbers of X-ray sources (filled squares) from Fig. [TJ 
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Fig. 9. — Computed T/N X b as a function of 7, showing scaling (see text). Computed curves 
for various values of T are closely bunched, as indicated. Overplotted are the positions of 
the globular clusters with significant numbers of X-ray sources (filled squares) from Fig. [H 
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Fig. 10. — Contours of constant N X b in the T — 7 (Verbunt parameters) plane. Overplotted 
are positions of Galactic globular clusters with significant numbers of X-ray sources detected 
in them (filled squares) from Fig. [TJ Nxb for each cluster is indicated by the number in 
parenthesis next to its marked position. 



